﻿using System;
using NGenerics.DataStructures.Mathematical;
namespace InvariantWM
{
    class DCT
    {
        public static double[,] directTransform(double[,] srcMatrix)
        {
            int n_rows = srcMatrix.GetLength(0);
            int n_cols = srcMatrix.GetLength(1);

            double[,] transformMatrix = calcDirectTransformMatrix(srcMatrix);

            //C = T*F*T_t;
            Matrix f = new Matrix(n_rows, n_cols, srcMatrix);
            Matrix t = new Matrix(n_rows, n_cols, transformMatrix);
            Matrix t_t = t.Transpose();
            Matrix c = (t.Multiply(f)).Multiply(t_t);
            
            return c.ToArray();
        }

        public static double[,] inverseTransform(double[,] srcMatrix)
        {
            int n_rows = srcMatrix.GetLength(0);
            int n_cols = srcMatrix.GetLength(1);

            double[,] transformMatrix = calcDirectTransformMatrix(srcMatrix);

            //F = T_t*C*T
            Matrix c = new Matrix(n_rows, n_cols, srcMatrix);
            Matrix t = new Matrix(n_rows, n_cols, transformMatrix);
            Matrix t_t = t.Transpose();
            Matrix f = (t_t.Multiply(c)).Multiply(t);

            return f.ToArray();
        }

        protected static double calculateBidimensionalInverseTransform(double[,] srcMatrix, int x, int y)
        {
            double coeff = 0.0;
            double alpha_u = 0.0;

            double alpha_v = 0.0;

            int n_rows = srcMatrix.GetLength(0);
            int n_cols = srcMatrix.GetLength(1);

            for (int u = 0; u < n_rows; u++)
            {
                for (int v = 0; v < n_cols; v++)
                {
                    if (u == 0)
                        alpha_u = Math.Sqrt(1.0 / n_rows);
                    else
                        alpha_u = Math.Sqrt(2.0 / n_rows);

                    if (v == 0)
                        alpha_v = Math.Sqrt(1.0 / n_cols);
                    else
                        alpha_v = Math.Sqrt(2.0 / n_cols);

                    double arg_u = ((2.0 * x) + 1.0) * Math.PI * u / (n_rows * 2);
                    double arg_v = ((2.0 * y) + 1.0) * Math.PI * v / (n_cols * 2);
                    double alpha = alpha_u * alpha_v;
                    coeff += alpha * srcMatrix[u, v] * Math.Cos(arg_u) * Math.Cos(arg_v);
                }
            }

            return coeff;
        }


        protected static double calculateBidimensionalTransform(double[,] srcMatrix, int u, int v)
        {
            double coeff = 0.0;
            double alpha_u = 0.0;
            double alpha_v = 0.0;

            int n_rows = srcMatrix.GetLength(0);
            int n_cols = srcMatrix.GetLength(1);

            for (int x = 0; x < n_rows; x++)
            {
                for (int y = 0; y < n_cols; y++)
                {
                    double arg_u = ((2.0 * x) + 1.0) * Math.PI * u / (n_rows * 2);
                    double arg_v = ((2.0 * y) + 1.0) * Math.PI * v / (n_cols * 2);
                    coeff += srcMatrix[x, y] * Math.Cos(arg_u) * Math.Cos(arg_v);
                }
            }

            if (u == 0)
                alpha_u = Math.Sqrt(1.0 / n_rows);
            else
                alpha_u = Math.Sqrt(2.0 / n_rows);

            if (v == 0)
                alpha_v = Math.Sqrt(1.0 / n_cols);
            else
                alpha_v = Math.Sqrt(2.0 / n_cols);

            return (alpha_u * alpha_v * coeff);



        }

        protected static double[,] calcDirectTransformMatrix(double[,] srcMatrix)
        {
            int n_rows = srcMatrix.GetLength(0);
            int n_cols = srcMatrix.GetLength(1);

            double[,] transformMatrix = new double[n_cols, n_rows];

            double alpha1 = Math.Sqrt(1.0 / n_cols);
            double alpha2 = Math.Sqrt(2.0 / n_cols);

            double cosine = 0.0;
            double arg = 0.0;


            //Compute first line
            for (int j = 0; j < n_rows; j++)
                transformMatrix[0, j] = alpha1;

            for (int i = 1; i < n_rows; i++)
            {
                for (int j = 0; j < n_cols; j++)
                {
                    arg = ((j + 0.5) * Math.PI * i) / n_rows;
                    cosine = Math.Cos(arg);
                    transformMatrix[i,j] = alpha2 * cosine;
                }
            }

            return transformMatrix;
        }

        protected static double calcDirectCoeff(int i, int j, int length, double alpha)
        {
            //double coeff = 0.0;
            double cosine = 0.0;
            double arg = 0.0;

            /*for (int j = 0; j < length; j++)
            {
                arg = ((j + 0.5) * Math.PI * i) / length; 
                cosine = Math.Cos(arg);
                coeff +=  cosine;                    
            }*/



            arg = ((j + 0.5) * Math.PI * i) / length;
            cosine = Math.Cos(arg);

            return (cosine*alpha);
        }

    }
}
